home *** CD-ROM | disk | FTP | other *** search
/ EnigmA Amiga Run 1998 July / EnigmA AMIGA RUN 29 (1998)(G.R. Edizioni)(IT)[!][issue 1998-07 & 08].iso / earcd / phase5 / ppcrelease / libm / tahoe / sqrt.s < prev    next >
Text File  |  1998-02-21  |  4KB  |  140 lines
  1. /*
  2.  * Copyright (c) 1987 Regents of the University of California.
  3.  * All rights reserved.
  4.  *
  5.  * Redistribution and use in source and binary forms, with or without
  6.  * modification, are permitted provided that the following conditions
  7.  * are met:
  8.  * 1. Redistributions of source code must retain the above copyright
  9.  *    notice, this list of conditions and the following disclaimer.
  10.  * 2. Redistributions in binary form must reproduce the above copyright
  11.  *    notice, this list of conditions and the following disclaimer in the
  12.  *    documentation and/or other materials provided with the distribution.
  13.  * 3. All advertising materials mentioning features or use of this software
  14.  *    must display the following acknowledgement:
  15.  *    This product includes software developed by the University of
  16.  *    California, Berkeley and its contributors.
  17.  * 4. Neither the name of the University nor the names of its contributors
  18.  *    may be used to endorse or promote products derived from this software
  19.  *    without specific prior written permission.
  20.  *
  21.  * THIS SOFTWARE IS PROVIDED BY THE REGENTS AND CONTRIBUTORS ``AS IS'' AND
  22.  * ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
  23.  * IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
  24.  * ARE DISCLAIMED.  IN NO EVENT SHALL THE REGENTS OR CONTRIBUTORS BE LIABLE
  25.  * FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL
  26.  * DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS
  27.  * OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION)
  28.  * HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT
  29.  * LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY
  30.  * OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF
  31.  * SUCH DAMAGE.
  32.  */
  33.     .data
  34.     .align    2
  35. _sccsid:
  36. .asciz    "@(#)sqrt.s    5.6    (ucb.elefunt)    10/9/90"
  37.  
  38. /*
  39.  * double sqrt(arg)   revised August 15,1982
  40.  * double arg;
  41.  * if(arg<0.0) { _errno = EDOM; return(<a reserved operand>); }
  42.  * if arg is a reserved operand it is returned as it is
  43.  * W. Kahan's magic square root
  44.  * Coded by Heidi Stettner and revised by Emile LeBlanc 8/18/82.
  45.  * Re-coded in tahoe assembly language by Z. Alex Liu 7/13/87.
  46.  *
  47.  * entry points:_d_sqrt        address of double arg is on the stack
  48.  *        _sqrt        double arg is on the stack
  49.  */
  50.     .text
  51.     .align    2
  52.     .globl    _sqrt
  53.     .globl    _d_sqrt
  54.     .globl    libm$dsqrt_r5
  55.     .set    EDOM,33
  56.  
  57. _d_sqrt:
  58.     .word    0x003c          # save r2-r5
  59.     movl    4(fp),r2
  60.     movl    (r2),r0
  61.     movl    4(r2),r1    # r0:r1 = x
  62.     brb      1f
  63. _sqrt:
  64.     .word    0x003c          # save r2-r5
  65.     movl    4(fp),r0
  66.     movl    8(fp),r1    # r0:r1 = x
  67. 1:    andl3    $0x7f800000,r0,r2    # r2 = biased exponent
  68.     bneq    2f
  69.     ret            # biased exponent is zero -> 0 or reserved op.
  70. /*
  71.  *                # internal procedure
  72.  *                # ENTRY POINT FOR cdabs and cdsqrt
  73.  */
  74. libm$dsqrt_r5:            # returns double square root scaled by 2^r6
  75.     .word    0x0000        # save nothing
  76. 2:    ldd    r0
  77.     std    r4
  78.     bleq    nonpos        # argument is not positive
  79.     andl3    $0xfffe0000,r4,r2
  80.     shar    $1,r2,r0
  81.     addl2    $0x203c0000,r0    # r0 has magic initial approximation
  82. /*
  83.  *                # Do two steps of Heron's rule
  84.  *                # ((arg/guess)+guess)/2 = better guess
  85.  */
  86.     ldf    r4
  87.     divf    r0
  88.     addf    r0
  89.     stf    r0
  90.     subl2    $0x800000,r0    # divide by two
  91.     ldf    r4
  92.     divf    r0
  93.     addf    r0
  94.     stf    r0
  95.     subl2    $0x800000,r0    # divide by two
  96. /*
  97.  *                # Scale argument and approximation
  98.  *                # to prevent over/underflow
  99.  */
  100.     andl3    $0x7f800000,r4,r1
  101.     subl2    $0x40800000,r1    # r1 contains scaling factor
  102.     subl2    r1,r4        # r4:r5 = n/s
  103.     movl    r0,r2
  104.     subl2    r1,r2        # r2 = a/s
  105. /*
  106.  *                # Cubic step
  107.  *                # b = a+2*a*(n-a*a)/(n+3*a*a) where
  108.  *                # b is better approximation, a is approximation
  109.  *                # and n is the original argument.
  110.  *                # s := scale factor.
  111.  */
  112.     clrl    r1        # r0:r1 = a
  113.     clrl    r3        # r2:r3 = a/s
  114.     ldd    r0        # acc = a
  115.     muld    r2        # acc = a*a/s
  116.     std    r2        # r2:r3 = a*a/s
  117.     negd            # acc = -a*a/s
  118.     addd    r4        # acc = n/s-a*a/s
  119.     std    r4        # r4:r5 = n/s-a*a/s
  120.     addl2    $0x1000000,r2    # r2:r3 = 4*a*a/s
  121.     ldd    r2        # acc = 4*a*a/s
  122.     addd    r4        # acc = n/s+3*a*a/s
  123.     std    r2        # r2:r3 = n/s+3*a*a/s
  124.     ldd    r0        # acc = a
  125.     muld    r4        # acc = a*n/s-a*a*a/s
  126.     divd    r2        # acc = a*(n-a*a)/(n+3*a*a)
  127.     std    r4        # r4:r5 = a*(n-a*a)/(n+3*a*a)
  128.     addl2    $0x800000,r4    # r4:r5 = 2*a*(n-a*a)/(n+3*a*a)
  129.     ldd    r4        # acc = 2*a*(n-a*a)/(n+3*a*a)
  130.     addd    r0        # acc = a+2*a*(n-a*a)/(n+3*a*a)
  131.     std    r0        # r0:r1 = a+2*a*(n-a*a)/(n+3*a*a)
  132.     ret            # rsb
  133. nonpos:
  134.     bneq    negarg
  135.     ret            # argument and root are zero
  136. negarg:
  137.     pushl    $EDOM
  138.     callf    $8,_infnan    # generate the reserved op fault
  139.     ret
  140.